ATLAS/ICESat-2 Monthly Gridded Sea Ice Thickness¶

This notebook highlights some simple code for visualizing monthly gridded sea ice thickness for three winter periods: winter 2018-19, 2019-20, and 2020-21.

import os
import xarray as xr
import numpy as np
import pandas as pd

# Remove warnings to improve display
import warnings 
warnings.filterwarnings('ignore')

# Import utils 
from utils.misc_utils import is2_interp2d, restrictRegionally, compute_winter_means
from utils.read_data_utils import read_is2_data, read_book_data
from utils.plotting_utils import plot_winter_means, plotArcticMaps, polarComparisonMaps, polar_fig_and_axes, interactive_lineplot, interactive_map

# Plotting
import hvplot.xarray
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
from textwrap import wrap
import matplotlib as mpl
import matplotlib.pyplot as plt

GeoAxes._pcolormesh_patched = Axes.pcolormesh
%config InlineBackend.figure_format = 'retina'
mpl.rcParams['figure.dpi'] = 200

1) Read data¶

Set desired date range¶

Our analysis looks at winter data, so we’ll only load data for the Northern Hemisphere’s winter season (Nov-Apr)

winter18_19 = pd.date_range(start="2018-11-01", end="2019-04-01", freq="MS")
winter19_20 = pd.date_range(start="2019-09-01", end="2020-04-01", freq="MS")
winter20_21 = pd.date_range(start="2020-09-01", end="2021-04-01", freq="MS") 
winter_months = winter18_19.append(winter19_20).append(winter20_21)
print(winter_months)

DatetimeIndex(['2018-11-01', '2018-12-01', '2019-01-01', '2019-02-01',
               '2019-03-01', '2019-04-01', '2019-09-01', '2019-10-01',
               '2019-11-01', '2019-12-01', '2020-01-01', '2020-02-01',
               '2020-03-01', '2020-04-01', '2020-09-01', '2020-10-01',
               '2020-11-01', '2020-12-01', '2021-01-01', '2021-02-01',
               '2021-03-01', '2021-04-01'],
              dtype='datetime64[ns]', freq=None)

Read in ICESat-2 data¶

If the data doesn’t already exist in the user’s local drive, the function read_is2_data will download the files from the google storage bucket.

ds = read_is2_data() # Read in data
ds = ds.sel(time=winter_months) # Get winter months
print(ds)

<xarray.Dataset>
Dimensions:            (time: 22, y: 448, x: 304)
Coordinates:
  * time               (time) datetime64[ns] 2018-11-01 ... 2021-04-01
    longitude          (y, x) float32 168.3 168.1 168.0 ... -10.36 -10.18 -9.999
    latitude           (y, x) float32 31.1 31.2 31.3 31.39 ... 34.68 34.58 34.47
    xgrid              (y, x) float32 -3.838e+06 -3.812e+06 ... 3.738e+06
    ygrid              (y, x) float32 5.838e+06 5.838e+06 ... -5.338e+06
Dimensions without coordinates: y, x
Data variables:
    projection         (time) float64 -2.147e+09 -2.147e+09 ... -2.147e+09
    ice_thickness      (time, y, x) float32 nan nan nan nan ... nan nan nan nan
    ice_thickness_unc  (time, y, x) float32 nan nan nan nan ... nan nan nan nan
    num_segments       (time, y, x) float32 nan nan nan nan ... nan nan nan nan
    mean_day_of_month  (time, y, x) float32 nan nan nan nan ... nan nan nan nan
    snow_depth         (time, y, x) float32 nan nan nan nan ... nan nan nan nan
    snow_density       (time, y, x) float32 nan nan nan nan ... nan nan nan nan
    freeboard          (time, y, x) float32 nan nan nan nan ... nan nan nan nan
    ice_type           (time, y, x) float32 nan nan nan nan ... nan nan nan nan
    ice_density        (time, y, x) float32 nan nan nan nan ... nan nan nan nan
Attributes:
    contact:      Alek Petty (alek.a.petty@nasa.gov)
    description:  Gridded Feb 2020 Arctic sea ice thickness and ancillary dat...
    reference:    Data doi: 10.5067/CV6JEXEE31HF. Original methodology descri...
    history:      Created 22/04/21

Set interpolation preferences¶

The dataset included uninterpolated data, or we can interpolate the data using linear interpolation or a simple nearest neighbor interpolation function. Because ICESat-2 doesn’t provide full monthly coverage, interpolating fills missing grid cells with a best guess based on surrounding data. This helps avoid sampling biases when performing time series analyses, with the cavaet that this interpolation method is subjective.

In order to define the interpolation bounds (so that we don’t try and interpolate over land, or other areas where there wouldn’t be sea ice!), we use the NOAA/NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration dataset to infer the location of the sea ice. By setting interpolation=True, the code calls an interpolation function that is defined in the misc_utils.py module.

%%time 

# Set interpolation preferences 
interpolate = True
interp_method = "nearest"

# Read book data. This contains the CDR data, which we'll use for interpolating 
book_ds = read_book_data()
book_ds = book_ds.sel(time=winter_months) # Get winter months

# Interpolate
if (interpolate == True): 
    print("Interpolating ICESat-2 data...")
    cdr_da = book_ds["cdr_seaice_conc_monthly"] # Get CDR data
    ds = is2_interp2d(ds, cdr_da, method=interp_method, interp_var='all')
    print("Complete!")
    
# Combine datasets 
ds = xr.merge([ds, book_ds])
ds = ds.drop_vars("projection")

Interpolating ICESat-2 data...

Complete!

---------------------------------------------------------------------------
MergeError                                Traceback (most recent call last)
<timed exec> in <module>

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/xarray/core/merge.py in merge(objects, compat, join, fill_value, combine_attrs)
    898         dict_like_objects.append(obj)
    899 
--> 900     merge_result = merge_core(
    901         dict_like_objects,
    902         compat,

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/xarray/core/merge.py in merge_core(objects, compat, join, combine_attrs, priority_arg, explicit_coords, indexes, fill_value)
    633 
    634     prioritized = _get_priority_vars_and_indexes(aligned, priority_arg, compat=compat)
--> 635     variables, out_indexes = merge_collected(
    636         collected, prioritized, compat=compat, combine_attrs=combine_attrs
    637     )

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/xarray/core/merge.py in merge_collected(grouped, prioritized, compat, combine_attrs)
    238                 variables = [variable for variable, _ in elements_list]
    239                 try:
--> 240                     merged_vars[name] = unique_variable(name, variables, compat)
    241                 except MergeError:
    242                     if compat != "minimal":

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/xarray/core/merge.py in unique_variable(name, variables, compat, equals)
    147 
    148     if not equals:
--> 149         raise MergeError(
    150             f"conflicting values for variable {name!r} on objects to be combined. "
    151             "You can skip this check by specifying compat='override'."

MergeError: conflicting values for variable 'ice_thickness_smoothed' on objects to be combined. You can skip this check by specifying compat='override'.

2) Restrict data to a given region¶

We’ve built a region mask for the Arctic into the dataset used for this book; see the data wrangling notebook for more information. The region mask is included as a coordinate in the dataset, allowing us to easily access the different regions. You can restrict the data by selecting keys corresponding to regions of interest.

We are often interested in a region called the Inner Arctic, which is defined as the combined area of the Central Arctic, Beaufort Sea, Chukchi Sea, E Siberian Sea and the Laptev Sea. Figure from Petty et al., (2020):

Inner Arctic region map

Each region is defined by a specific integer key value. We can select the region/s defined by the key/s by using the restrictRegionally function, defined in the misc_utils module. Here, we’ll restrict the data to the Inner Arctic by selected the key values [10,11,12,13,15], corresponding to the five regions described above.

regions = pd.DataFrame({"labels":ds.region_mask.attrs["labels"]}, index=ds.region_mask.attrs["keys"])
print(regions)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/var/folders/8v/cr06mz0n3bjd5jr12_9v6n200000gn/T/ipykernel_81929/1017187880.py in <module>
----> 1 regions = pd.DataFrame({"labels":ds.region_mask.attrs["labels"]}, index=ds.region_mask.attrs["keys"])
      2 print(regions)

~/opt/anaconda3/envs/icesat2_book/lib/python3.8/site-packages/xarray/core/common.py in __getattr__(self, name)
    237                 with suppress(KeyError):
    238                     return source[name]
--> 239         raise AttributeError(
    240             "{!r} object has no attribute {!r}".format(type(self).__name__, name)
    241         )

AttributeError: 'Dataset' object has no attribute 'region_mask'

innerArcticKeys = [10,11,12,13,15] #Inner Arctic
ds_innerArctic = restrictRegionally(ds, innerArcticKeys)

Regions selected: Inner Arctic

On the map below, you can see how the data has been cut to show only gridcells in the Inner Arctic!

date = "Feb 2019"
fig, axes = polar_fig_and_axes(hemisphere="nh", num_rows=1, num_columns=2, lat=55, show_map_features=True)
im1 = ds["ice_thickness"].sel(time=date)[0].plot.pcolormesh(ax=axes[0], x='longitude', y='latitude', vmin=0, vmax=4, extend='both', transform=ccrs.PlateCarree(), add_colorbar=False)
im2 = ds_innerArctic["ice_thickness"].sel(time=date)[0].plot.pcolormesh(ax=axes[1], x='longitude', y='latitude', vmin=0, vmax=4, extend='both', transform=ccrs.PlateCarree(), add_colorbar=False)
axes[0].set_title("Arctic sea ice thickness ("+date+")")
axes[1].set_title("Inner Arctic sea ice thickness ("+date+")")
plt.colorbar(im1, cax=fig.add_axes([0.93, 0.326, 0.025, 0.352]), extend='both', label="sea ice thickness (m)")
plt.show()

3) Plot mean sea ice thickness over time¶

We’ll use the function plot_winter_means defined in the plotting_utils module to generate a line plot of mean sea ice thickness by month. This function also allows you to plot the uncertainty as shading around the line. Sea ice thickness is provided as a data variable in the ATLAS/ICESat-2 Monthly Gridded Sea Ice Freeboard dataset.

# Do you want to use the Inner Arctic restricted data or the entire, unrestricted region? 
use_innerArctic = True

# Grab sea ice thickness and uncertainty from the dataset 
if use_innerArctic == True: 
    thickness = ds_innerArctic["ice_thickness"] # Sea ice thickness
    thickness_unc = ds_innerArctic["ice_thickness_unc"] # Sea ice thickness uncertainty 
else: 
    thickness = ds["ice_thickness"] 
    thickness_unc = ds_innerArctic["ice_thickness_unc"]
    
# Descriptive info regarding whether or not the data has been interpolated
if interpolate==True: 
    interp_descrip = " (interpolated)" # Description to add to plot titles 
else: 
    interp_descrip = "" # empty string

Compute mean monthly winter sea ice thickness¶

ICESat2 has varying coverage by month based on its oribal path and cloud coverage. In order to avoid the mean being artificially inflated or deflated due to varying coverage over the course of the season, we ensure that we’re only computing the mean in regions that have good coverage. This is looped into the compute_means function; see the function in the misc_utils module for more information.

mean_winter_thickness = compute_winter_means(thickness)
display(mean_winter_thickness)

time period	Winter 2018-19	Winter 2019-20	Winter 2020-21
month
Sep	NaN	1.353715	1.161614
Oct	NaN	1.270545	1.052572
Nov	1.206516	1.202967	0.927755
Dec	1.429245	1.327576	1.131299
Jan	1.647857	1.615521	1.456750
Feb	1.890510	1.848658	1.697177
Mar	2.021039	1.952370	1.870251
Apr	2.132214	1.994546	1.979356

Static plot¶

Here, we’ll generate a lineplot to show the winter means and the associated uncertainty in the measurement.

##### STILL WORKING ON UPDATING THIS FUNCTION TO REFLECT THE BUSINESS OF NANNING GRIDCELLS ACROSS MONTHS. 
##### WOULD BE GOOD TO COMBINE THE PLOTTING FUNCTIONALITY WITH THE compute_winter_means FUNCTION IN THE misc_utils MODULE

#start_year = str(winter_months[0].year) 
#end_year = str(winter_months[-1].year)
#plot_winter_means(da=thickness, start_year=start_year, end_year=end_year, 
#                  show_uncertainty=False, 
#                  title="Mean monthly sea ice thickness"+interp_descrip) 

Interactive lineplot¶

This function uses the package hvplot to generate an interactive lineplot, allowing the user to see the values of each point on the line

interactive_lineplot(mean_winter_thickness, ylabel="mean sea ice thickness (m)", title="ICESat2 mean monthly sea ice thickness")

4) Map mean winter ice thickness¶

Using the cartopy mapping package for python, we overlay mean ice thickness for each winter season on a basemap of the Arctic.

NOTE: ICESat-2 only started collecting data in November 2018. The periods for the mean has been shortened to Nov-Apr, and the line plots for winter 2018-19 will show less data, since no data was available for September or October 2018.

Compute gridcell mean sea ice thickness¶

Here, we’ll compute the mean value in each gridcell for the period of Nov-Apr for each winter season.

# Compute mean and add descriptive attributes
winter1_mean = thickness.sel(time=winter18_19).mean(dim="time").assign_coords({"time":"Nov 2018 - Apr 2019"})
winter2_mean = thickness.sel(time=winter19_20[2:]).mean(dim="time").assign_coords({"time":"Nov 2019 - Apr 2020"})
winter3_mean = thickness.sel(time=winter20_21[2:]).mean(dim="time").assign_coords({"time":"Nov 2020 - Apr 2021"})

Plot winter means on a map of the arctic¶

Here, we’ll use the function plotArcticMaps defined in the plotting_utils module to generate the maps.

# Generate maps of each mean    
means_all = xr.concat([winter1_mean, winter2_mean, winter3_mean], dim="time")# All means in one dataset
plotArcticMaps(means_all, minval=0, maxval=4, title="Gridcell mean sea ice thickness"+interp_descrip)

Compare the means for each winter¶

The rightmost plot shows gridcell difference in the means, allowing you to better understand spatial patterns in differences between two winter seasons. To generate these plots, we use a function named polarComparisonMaps, which takes two datasets as its input. See the plotting_utils module to see the function and more details on its utility and construction. Here, we’ll feed the function the gridcell means for two winter seasons, and some descriptive labels and such in the function arguments to improve the figure output.

y="latitude"
x="longitude"
colorbar=True
clabel="sea ice thickness (m)"
vmin=None
vmax=None
cmap="viridis"
title=None
hemisphere="nh"
lat=60
max_lat=90

pl1 = winter1_mean.hvplot.quadmesh(y=y, x=x, 
                   colorbar=False, clim=(0,4), cmap=cmap, 
                   title=winter1_mean.time.values.item(), features=["coastline"],
                   projection=ccrs.NorthPolarStereo(central_longitude=-45), geo=True, project=True, 
                   frame_width=215, ylim=(lat,max_lat), xlim=(-180,180), 
                   dynamic=False)
pl2 = winter2_mean.hvplot.quadmesh(y=y, x=x, 
                   colorbar=False, clim=(0,4), cmap=cmap, 
                   title=winter2_mean.time.values.item(), features=["coastline"],
                   projection=ccrs.NorthPolarStereo(central_longitude=-45), geo=True, project=True, 
                   frame_width=215, ylim=(lat,max_lat), xlim=(-180,180), 
                   dynamic=False)
pl3 = winter3_mean.hvplot.quadmesh(y=y, x=x, 
                   colorbar=True, clim=(0,4), clabel=clabel, cmap=cmap, 
                   title=winter3_mean.time.values.item(), features=["coastline"], 
                   projection=ccrs.NorthPolarStereo(central_longitude=-45), geo=True, project=True, 
                   frame_width=215, ylim=(lat,max_lat), xlim=(-180,180), 
                   dynamic=False) 

pl1

import holoviews as hv

y="latitude"
x="longitude"
colorbar=True
clabel="sea ice thickness (m)"
vmin=None
vmax=None
cmap="viridis"
title=None
hemisphere="nh"
lat=60
max_lat=90

pl1 = winter2_mean.hvplot.quadmesh(y=y, x=x, 
                   colorbar=False, clim=(0,4), cmap=cmap, 
                   title=winter2_mean.time.values.item()+" mean", features=["coastline"],
                   projection=ccrs.NorthPolarStereo(central_longitude=-45), geo=True, project=True, 
                   frame_width=215, ylim=(lat,max_lat), xlim=(-180,180), 
                   dynamic=False)
pl2 = winter3_mean.hvplot.quadmesh(y=y, x=x, 
                   colorbar=True, clim=(0,4), clabel=clabel, cmap=cmap, 
                   title=winter3_mean.time.values.item()+" mean", features=["coastline"],
                   projection=ccrs.NorthPolarStereo(central_longitude=-45), geo=True, project=True, 
                   frame_width=215, ylim=(lat,max_lat), xlim=(-180,180), 
                   dynamic=False)
diff = winter3_mean - winter2_mean
diff.attrs["long_name"] = "sea ice thickness gridcell difference"
pl3 = diff.hvplot.quadmesh(y=y, x=x, 
                   colorbar=True, clim=(-1.5,1.5), clabel="gridcell difference", cmap="coolwarm", 
                   title="winter 2020-21 - winter 2019-20", features=["coastline"], 
                   projection=ccrs.NorthPolarStereo(central_longitude=-45), geo=True, project=True, 
                   frame_width=215, ylim=(lat,max_lat), xlim=(-180,180), 
                   dynamic=False) 


data_pl = (pl1+pl2+pl3).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
display(data_pl)
hvplot.save(data_pl, 'mean_diff.html')

# Winter 2018-2019 and Winter 2019-2020 differences 
polarComparisonMaps(data1=winter1_mean, label1="Winter 2018-19", 
                    data2=winter2_mean, label2="Winter 2019-20", 
                    vmin=0, vmax=4, vmin_diff=-1.5, vmax_diff=1.5, 
                    title="Mean winter sea ice thickness (November -> April)")

# Winter 2019-2020 and Winter 2020-2021 differences 
polarComparisonMaps(data1=winter2_mean, label1="Winter 2019-20", 
                    data2=winter3_mean, label2="Winter 2020-21", 
                    vmin=0, vmax=4, vmin_diff=-1.5, vmax_diff=1.5, 
                    title="Mean winter sea ice thickness (November -> April)")

# Winter 2018-2019 and Winter 2020-2021 differences 
polarComparisonMaps(data1=winter1_mean, label1="Winter 2018-19", 
                    data2=winter3_mean, label2="Winter 2020-21", 
                    vmin=0, vmax=4, vmin_diff=-1.5, vmax_diff=1.5, 
                    title="Mean winter sea ice thickness (November -> April)")

Generate interactive maps¶

This allows you to scroll over individual gridcells and zoom in on regions of interest

import holoviews as hv
import hvplot.xarray

from utils.plotting_utils import compute_vmin_vmax
y="latitude"
x="longitude"
colorbar=True
vmin=None
vmax=None
clabel=""
cmap="viridis"
title=None
hemisphere="nh"
lat=60
width=600
height=350 

da = thickness.sel(time=winter18_19)
vmin=0
vmax=4
clabel="sea ice thickness (m)"

if hemisphere=="nh": 
    max_lat = 90
    projection = ccrs.NorthPolarStereo(central_longitude=-45)
elif hemisphere=="sh": 
    max_lat = -90 
    projection = ccrs.SouthPolarStereo()
    if lat > 0: 
        lat = -lat
    
vmin_auto, vmax_auto = compute_vmin_vmax(da)
clim = [vmin,vmax]
if vmin is None: 
    clim[0] = vmin_auto
if vmax is None: 
    clim[1] = vmax_auto 

help(da.hvplot.quadmesh)

da.hvplot.quadmesh(y=y, x=x, 
                   colorbar=colorbar, clim=tuple(clim), clabel=clabel, cmap=cmap, 
                   title=title, features=["coastline"],
                   projection=projection, geo=True, project=True, 
                   ylim=(lat,max_lat), responsive=True,
                   dynamic=False)

interactive_map(thickness.sel(time=winter19_20), vmin=0, vmax=4, clabel="sea ice thickness (m)")

interactive_map(thickness.sel(time=winter20_21), vmin=0, vmax=4, clabel="sea ice thickness (m)")

rasm = xr.tutorial.open_dataset('rasm').load()



rasm.hvplot.quadmesh(
    'xc', 'yc', crs=ccrs.PlateCarree(), projection=ccrs.PlateCarree(),
    ylim=(0, 90), cmap='viridis', project=True, geo=True,
    rasterize=True, coastline=True, frame_width=800, dynamic=False,
)

ICESat-2 Sea Ice Thickness Data Analysis and Visualization